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Abstract 

In this paper, we focus on the analytical expressions of both three-point and four-point integrals 
for the case of a small Gram determinant. We also investigate the numerical efficiency of n- 
point integrals for n > 5. One-loop five-point and higher-point integrals are crucial for accurate 
predictions of LHC physics, such as Higgs searches involving background process with three or more 
final states. The expressions of five-point and six-point integrals have been calculated in Ref. [l| 
using the fact that the loop momentum can be expressed as a linear combination of the external 
momenta. We first present a proof of the results of Ref. [l| using a new reduction formula. NLO 
calculations typically suffer numerical instabilities in specific phase space regions when individual 
kinematic variables are small; this leads to a vanishing Gram determinant. In these problematic 
regions, we introduce a technique that resolves the numerical instability analytically instead of 
using a numerical iteration. Explicit expressions of three-point and four-point integrals in the 
limit of small Gram determinants are provided by this new method. Furthermore, we demonstrate 
that one-loop n-point integral (with n > 6) can always be reduced to six (n— l)-point integrals in 
general. This dramatically reduces the CPU time for numerical computations compared to other 
methods. Additionally, the numerical uncertainty originating from the computation of the higher 
dimensional Cayley matrix can be removed. We present general reduction formulas for five-point 
and higher-point scalar, vector, and tensor integrals at the one-loop level. 
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FIG. 1: General one- loop n-point diagram for the standard form of T™ ... with pij = Pi — Pj- The 
arrows show the direction of momentum rather than the particle's current. 

I. INTRODUCTION 

Recent experiments have produced precise predictions for Standard Model (SM) observ- 
ables in both the electroweak (EW) and QCD sectors. Many of these observables involve 
multi-particle final-states. To make accurate NLO predictions for these observable, it is nec- 
essary to compute loop diagrams with many legs. The key ingredient of these calculations 
is the n-point loop integral. Accurate theoretical predictions therefore rely on accurate and 
efficient evaluation of the various n-point integrals which we encounter. 

In this paper, we study one-loop integrals that are essential for the NLO radiative correc- 
tions. We are particularly interested in one-loop integrals for processes with three or more 
final-state particles. The various n-point integrals necessary have been calculated partially 
or completely j^rB]. If the conventional Passarino-Veltman reduction Q scheme is used, 
we must address the difficulty of a vanishing Gram determinant. Also, if we are comput- 
ing heavy quarks, our results must be solved for general mass values, and this drastically 
complicates the kinematics. 

This paper is organized as follows: In Sec. [Til we prove the result of the Ref. [l| for n- 
point integrals (n > 5) with algebraic formalism. In Sec. lIII[ we drive the explicit expressions 
of three-point and four-point integrals for the case of a vanishing Gram determinant. In 
Sec. IIV[ the expressions of n-point integrals with n > 5 will be provided. In Sec. IIV[ the 
numerical improvement will be summarized. 

II. FORMALISM 

A. For the case of linear dependence on external momenta 

We begin with some standard definitions for the one-loop scalar n-point integral following 
Ref. The n-point scalar integral in standard form is: 




(1) 



2 



with the denominator factors: 

N k = (q + p k ) 2 ~ml + ie , k = 0, • • • , n - 1 , p = . (2) 

In general, we also use alphabetic notation for the n-point functions; T 1 = A, T 2 = B, 
T 3 = C, T 4 = D, and so on. For the reduction formula, the integral with one denominator 
factor omitted is very useful; (k) denotes that k th denominator N k is omitted in Eq. (JT]): 

^(*h^/^=-^4 — =- • (3) 



m 2 J ^N ---N k _ 1 N k+1 ---N n _ 1 

As shown in Ref p], it is assumed that Tq can be spanned by T^ _1 (A;), which we will prove 
with a different method: 

n— 1 

i3 , = X) a * i r 1 (0 • ( 4 ) 

8=0 

If we find all a, coefficients, then we can reduce any n-point scalar integral to a number of 
(n — l)-point scalar integrals. 

Let us start from the integrand of Eq. (pQ) with unknown dj coefficients: 

+ a r Ar at r • • • + —7 — — — — • [0) 



N ---N n „t Nf..N n -i N N 2 -..N n -i N ---N n _ 2 N Q - • ■ N n -\ 

This must be true for arbitrary q^. Comparison of numerators in both sides gives: 

n— 1 n— 1 



1 = ^ a^iV; = ^ a* (q 2 - 2g u pf + ^ 2 - m 2 ) 

i=0 i=0 

\ n-1 

53 ai I - 2g M «" + (Pi ~ m ?) a * > ( 6 ) 



q 2 



where k^ 1 is: 



, i=0 / i=Q 



n-1 



i=0 



In order for Eq. ([6]) to always be true for all q^, the coefficients of q^ and q 2 must vanish, 
while the constant term should be equal to 1. This implies: 



n-1 



5> = > ( § ) 



i=0 

n-1 



5^-77^ = 1 , (9) 



i=0 



= . (10) 
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Most importantly, we should find a prerequisite for the validity of reduction formula assumed 
in Eq. (j3J) by taking a careful look at the condition of Eq. ( TTUj) : 

k m = p^ax + p%a 2 H h p^an-! = . (11) 

If a set of vectors, {pf,i — 1, • • • , n — 1}, are linearly independent, then every coefficients 
of those vectors, ao, ai, • • • , a n _i, are all zero, which is a trivial solution. In order to have 
non-trivial solution, the set of those vectors should, therefore, be linearly dependent, where 
the total number of vectors should be greater than the dimension of vector space D. Con- 
sequently, n — 1 > D must be satisfied. In four-dimensional vector space, n > 5. (This 
constraint is not valid for the linearly independent case in Sec III Bp . 

Next, we should consider whether a system of linear equations, Eq. ([8]), Eq. (J9j, and 
Eq. ffTUj) . is solvable or not. Eq. ([8]) is equivalent to one of Eq. (fTOl) because Eq. ffTOl) doesn't 
change when we apply Eq. ([8]) to Eq. ffTUl) . Therefore, the total number of equations is five 
while the number of unknowns is n. Consequently n > 5 for the underdetmined system. 
This means that the reduction formula, Eq. (jl]) is valid only for the more than five point 
scalar integral. 



/t M = in Eq. ffTUj) can be modified as a scalar equation by dotting 



n-1 



2p E ■ k = 2^2(p E -pi)ai = 



(12) 



i=0 



where £ can be one of {0, 1, • • • , n — 1}, and a factor 2 is introduced for later use. 
In an attempt to derive a concise matrix form, let us begin with Eq. f|T2|) : 



n-1 



n-1 



= 2 J2(Pt ■ Pi) a i = 2 Pi a £ + 2 (P- 



rVi) a i 



t=i 



n-1 



n-1 



pja e + ( 1 + mja e - ^ (P? ~ m ?) a * I + 2 $Z ^ rPi 

i=0,i^e ) i=0,i^£ 



n-1 

1 + 2mja£ — 



(P 2 i ~ m i + Pi ~ m \ ~ 2 PrPi) ^ 



n-1 



Pi) 



i=0 



where is defined as pj — pj. To reach the second line from the first, we applied Eq. 
We also used a i = 0- 

Therefore, by varying i from to n — 1: 



n-1 



n-1 



3=0 j=0 

Here, the matrix Y is called Cayley matrix. Now a system of equations becomes: 
/ 

O O Q O O O 

ai 



(13) 



m\ + vn\- pl 



™o + m n-i ~:Pon-i\ ( a \ 
m\ + m 2 a _ x -p\ n _ x 



WLi + m l - Pl-i o 



2 K-x 



J 



1 



W-i/ 
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FIG. 2: Diagrammatic description of the reduction (Eq. (|24|) ) introduced in Sec. Ill Bl (n + l)-point 
integral reduces to n-point integral. pg is chosen to be p\, and p n = p n o should be (anti)parallel 
to pg. In order to avoid unphysical singularities, external legs pio, p n o and the first and last loop 
propagators should share the same mass term. 

If a determinant, \Y\, does not vanish, the above matrix equation can be completely solved 
for the aj coefficients. Consequently, the n-point scalar integral satisfies the following reduc- 
tion formula: 

n-1 

T " = ^a 4 T - 1 W, n>5 , (14) 

i=0 

with 

n-1 

Oi = -X)( y "%> (< =0,..-,n-l) . (15) 

3=0 

This result is exactly identical to that of Ref. [l[ for n-point scalar integrals with n > 6. 

For discussion, let us consider 6-point scalar integral. The reduction formula is based on 
solving a system of equations including k m = 0. 

k m = aip^ + a 2 p^ + a 3 p% + a^l + a 5 p£ = . 

Physically and mathematically, it is always true that any five vectors are linearly dependent 
in a four-dimensional vector space if all of five momenta does not vanish. Difficulties would 
appear when, for instance, all components of p± are so small to be negligible, which may 
correspond to soft gluon radiation from a hexagon diagram; in this specific case, p%, ■ ■ ■ ,p§ 
could be linearly independent so that reduction formula, Eq. (TT4"j) is not safe. The remedy 
of which will be provided in Sec. Ill Bt (The resultant expression is identical). 

B. For the case of linear independence on external momenta 

The linear independence on p^s prohibits the use of compact formula that is derived at 
previous section. However, we can introduce p n as a n th denominator N n to Eq. §5§ in order 
to break down the independence of Pi, • • • ,p n -i- 

Pn = xpi . 
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With N n , we can apply the same identity as described in Sec. Ill Al 

N ---N n ^N n N ---N n ^N n ' y ' 
where aj = — ^j =0 O^ -1 )^-- Multiplying q 2 — in both sides gives: 



N ---N n ^N n N ---N n ^N n 



Vg^ . (16) 



While x = gives a vanishing Cayley determinant for the (n + l)-point loop, the limiting 
case of x — > reduces this (n + l)-point into n-point integral as shown in Fig. |2j For small 
x, the unphysical external leg of po n (= xpoi) can be absorbed by the external leg of pio 
when mass term is shared. Besides, a loop propagator of q + xp± can also be absorbed by 
the external leg of q with the shared mass term. The fact that the fermion current is always 
involved in a loop guarantees the existence of a common mass term. The identity Eq. ( II 6p 
supposes that a fermion current consists of p nn -i, q and p w . The mass of the particle on 
an external leg of po n is assumed to be m , which is also the mass of n th propagator N n as 
illustrated in Fig. Employing the limit x — > takes it back to the result of the linearly 
dependent case in Sec. Ill Al 

r» = lim I q2 - ml = lim I ^oV-^K^ 
<*+°Jd N ---N n ^N n x^J D N ---N n ^N n 

nT'.n = E*r n - 1 (o • ( 17 ) 



i=0 



Hence, the previous reduction formula described in Sec. Ill Al is now extended to the region 
where external momenta are linearly independent. 

In order to verify the reduction procedure illustrated in Fig. [5] numerically, we should 
expand to second-order in x. For this goal, let us start from vector integrals: 



7i — 1 n—1 



(19) 



T n + 2b e J2Pi^T[ l = J2 a ^ Tn ~ 1 ^ ■ ( 18 ) 

This reduction can be easily derived by considering the equation as follows: 

2b t q-p t + 1 YaZq ai N i 
Nq - ■ ■ N n _i ~ No - ■ ■ N n -i 

Coefficients a^'s and be must be solutions of a system of equations shown below: 

n—1 n—1 

^2 a. = , ^2(P 2 ~ ml)ai = 1 , a x p^ H a n ^ n _ x = b e p e . 

i=0 i=0 

The fact that the p^s are linearly independent implies that all a^s vanish except and bf. 

a + bi = , — mla + (p 2 — m 2 )ai = 1 , 
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FIG. 3: Comparison between Eq. (|14p and Eq. (|24p for five-point scalar integral (Eq) at different 
energy scales (y/s). Those two reductions generate identical numerical results which show the 
reduction formula described in Sec. Ill Al can be extended to the region where momenta are linearly 
independent. The Feynman diagram for this computation is illustrated at Fig. HI 



which gives: 



b £ = ai = -a = -jr • (20) 

h 



With these coefficients, and Eq. ( 1191) . we can derive a reduction formula for the vector 
integral by varying I — 1 , • • • , n — 1 : 

n-l 

i? = J2 [t^O) - Tn_1 (°) - > ( 21 ) 

where the Gram matrix, Z, is defined as = 2pi-pj, and = -m-+mg. This reduction 
formula is exactly the same form of the conventional Parssarino-Veltman reduction. 
For later use, the identity, Eq. (fT9|) . leads to: 



Nq - ■ ■ N n _i N ---N n ^ N - ■ ■ N n _i Nq - ■ ■ N n _i 

Now, we are ready to return to the scalar integral. By means of multiplying Eq. (JSJ) by 
(q 2 — ttIq) after the shift of n — >■ n + 1, the identity for the T n+1 scalar integral is: 

This is the same form of the identity of Eq. f fLSj) . Multiplying both sides with N n / (q 2 — mj) 
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gives: 



E' %/Vj QqA^q a w (g 2 + x 2 pj + 2xq-p £ - rog) 

1=1 

E' ajNj (oq + q n )A^ a w (x 2 p|) xa^q-pi) 

N ---N n ^ N ---N n ^ W ---W n _! iV Q ---W n _ a ' 

Taking integration with the limit a; — >■ 0, and the help of Eq. (}22l) . leads to the reduction 
formula: 

n-l 

T n = ^ b i T n ^ 1 {i) , (24) 

i=0 

where the frj's are: 

_ ^ Qj + fea;a n + & (1 - x)a n 
x^o 1 - a n (x 2 pj - xf e ) 

Numerical computation of the above reduction formula is identical to the general reduction 
formula of Eq. ( TT4"|) as shown in Fig. [3J 



C. For the one- loop vector and tensor n-point integrals 

Reduction for the vector and tensor integrals can be determined without difficulty. Al- 
though we multiply Eq. dSJ) and Eq. (fl6l) with q^ 1 ■ ■ ■ q^ r on both sides, the identity is still 
valid for all q^. Hence, the following formula can be obtained: 

T n _ (27T/i) 4 - D f D q^---q„ r 



o- 2 y - a jv iVi...jv n _i 

n-l 



E^irW • (25) 



i=0 



s 



Tensor integrals can be decomposed to scalar coefficient functions: 



[r/2] n -l 

(26) 



fe=0 11,12, ••• ,lr=l 



2 A- 



where we used combitoric notation < p *i'" p i r K which denotes that; First, pick out 2k out 

of {pf 1 • • -pf r }. Second, switch picked-out p^p" • • • to g^ v g" ■ ■ ■ as a sum of all possible 
combinations of Lorenz indices followed by leftovers of p's. Finally, the sub-indices of 
momentum i must be kept same order starting from i\ to i r -2k which is just dummy indices 
for summation, [r/2] is the maximal integer lass than or equal to r/2. Note that we 
apply possible combinations only for Lorenz indices, not for the momentum subscripts. As 
examples, two-point functions are decomposed: 



P v P 



0*X + 9 W VI + <f% 



which gives: 



Bp = PifiBi, 

B^v = V\pV\vB\\ + g^ v BQQ, 

B^vp = V\nV\vV\pB\\\ + (g^Pip + QupVxu + gv P V\p) -Booi, 
Bpupa = v\iiV\vV\pV\aB\w\ + {gp P g p( j + gu P g^a + gp^gua) B 0000 

+ (gpuPipPia + gppPiaPiu + gpuPiuPip + gupPiuPip + gpaPivPip + swpipPi^) -Boon- 

(27) 

Using the decomposition of Eq. ( 126]) . (T n ^ 1 ) fll '" tlr (k) that is defined in Eq. ([3]) becomes: 

fc=0 ll,— ,«r = l 2k 

[r/2] n-l 

.(0, (28) 



L ' J (ill jLtr^l 

[T^r^ (o)=e e r i 2ifc <p } :z s*- 

fc = it,— ,i r = l 2k 



■ % r-2k 



with, 



0*3 



i when z < j 
i — 1 when z > j ' 

5 fj = 1 - 5 tj . (29) 



We can directly derive the vector and tensor coefficient functions by comparing Eq. ( |26|) 
with Eq. (I2"8"|) since they are connected by the reduction formula of Eq. ([25]) : 



n-l 



^00—00 ii— i r _2! ~~ a 0^oo---OOx'i---i r _2i (^) ~^ ^ ] a k$iik ' " " ^iiA;^no...oo x'i — i r -2! (^^ 



2? 



fe=l 
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FIG. 5: Numerical safety of three-point functions : C..\p\,p\,t, m\, m\, m|]=C...[0, 0, t, 0, Myy, M§] 
with M w = 80.425 GeV, M z = 91.1876 GeV, And t varies from -0.0001 to 0. The conven- 
tional Passarino-Veltman reduction shows numerical instability in small kinematic region (red 
curve, computed via LoopTools [9J), while blue curve which generated from the method in Sec. IIIII 
(Eq. (|43|) . (|46|) ) is numerically stable. The Feynman diagram for this computation is shown in 
Fig.! 



III. THREE POINT AND FOUR POINT FUNCTIONS, C, D. FOR THE RE- 
GION OF SMALL KINEMATIC VARIABLES 

In Sec. [Til we studied the algebraic formalism for the reduction of one-loop integrals. 
The coefficient functions of four-point or less-point integrals reside in inner-product space 
where the dimension of vector space is four in general. Hence, the number of momentum 
involved in the loop is always less than four. This mostly guarantees the dependence of 
momenta. Expressions for one- (A), two- {B), three- (C), and four-point (D) functions 
have already been known in an analytical form 0, However, we need the solutions to 



fix the problem that appear in the numerical computation, such as the case of the vanishing 
Gram determinant. In Ref. Q, the technique of numerical iteration is suggested for the 
small Gram determinant. 

The method in Sec. [Ill provides the analytical approach for the case of small kinematic 
variables. If Pi-pj = and pf = for all i and j, then the identity of Eq. ([5]) can be greatly 
simplified. Furthermore, every n-point functions can be reduced to a sum of two-point 
functions: 

1 ^ A 



1 = 1 
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After some simple algebra, the Ai : s are easily determined: 



n-l 



A = nB -w \ iT , X> = • ( 31 ) 

{ml -mfj ^ 



If the Pi-p/s are small enough, the loop propagator can be expanded around zero using 
the small variables, z = p 2 , z = 2q-p and z = z + z: 



z + z (z + z) 2 



(q + p) 2 — m 2 q 2 — m 2 (q 2 — m 2 ) 2 (q 2 — m 2 ) 3 

= E (^rpysr = E = E h (-^)* < Ar ') • (32) 

fc=0 v y fc=0 fc=0 

where the denominator M = q 2 — m 2 and the differential operator d = 8/ dm 2 . Therefore, 
with the small values of = p 2 , zi = 2q-p i: Zi = z% + 5», the identity of Eq. (jSJ) can be 
expanded as follows: 



1 



N --- iV n _i 



1 

2| -~ j 3 ; 



1 - Zidi + —^ZiZjdidj - —.ZiZjZkdidjdk + 



1 



M • ■ ■ M n _x 



^A^Vi ~~ ^""^ ^ j + ^'''iVoW + 2! ^ fc + 2 ^ fc + ~ Zj ~ Zk) ^ dkAi "> n~n: 
- 7y {zjZ k z e + 3zjZ k z e + 3zjZ k ze + Zjhze) (djd k d e Ai) H 

2° f^j A ~ % ( d j>Ai) + T^ZjZk (djd k Ai) - ^ZjZ k z e (djd k d t Ai) H J — 

/ 1 2 3 \ 
+ 2 1 f - — (djAi) + — 4 (djdkAi) - gfkZi (djdkdeAi) H j p^ 



N Ni 



2 2 (9jd k A) - ^z t (djdkdeAi) p 



... | „M v%&L_ J 

' ' l PjPk N N { 



~ *' (0) iVoiV, + A >« iVoiV, + A '( 2 ) iVoiV, + " ' ' l ^ j 

where repeated sub- indices imply the summation from 1 to n — 1, and coefficient functions 
J 7 are defined implicitly 

By integrating Eq. (133|) . the n-point scalar integral can be expressed in terms of two-point 
functions: 

r o " = T ij( o)^ + ^ (1) 4 + ^ 2) ^ + ... > (34) 
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where, 



Bi = B...(plm 2 ,ml 



(35) 



Vector and second rank tensor coefficient functions, T™, can be obtained from Eq. (I3"3"|) by 
multiplying or q^q u : 



Tm nr~ rfi i rtr-v r>i i nr^P r>i i 
„ = J-i,(o) B ii + Am^v + A,(2) i Vp + " ' 

L pv - ^,(0)-D MI , + -Ti^ixvp + ^i,[2) U fJ,vpT + 



(36) 



When two masses are identical (say, m a = m&), it is necessary to take an additional 
derivative with respect to m\ in Eq. (|33|) : 



1 



M ---M a 2 ---M n _ 2 



-d a 



1 



Mo • • • M a ■ ■ ■ M n _ 2 



(37) 



Similarly, if n masses are the same, we should take n-th order of derivative for the mass. 

In Appendix |AJ we present the explicit expressions of the 5-functions and their deriva- 
tives. 

In Eq. fl34l) . n-point functions are expressed by two-point (B) functions. Alternatively, 
we can expand Eq. fl33l) in terms of three-point (C) functions when two masses are identical. 
As an example, four-point integrals with m 2 = m% for the case of a small Gram determinants 
yields: 



Nn - ■ ■ 



1 
1 

A3 



1 



N Q ---N 2 

% 1 



-i,(0) 



— 1— Zp» 

+ ^i, 1 



■fMU 



i,(0) 



iVoiV, ' ' ^AW ^MjVoJVi 



+ 



A A,A 3 



+ 



Therefore, 



where, 



L)q — I -J- i,(0)^o + ^ + ^ 1,(2)°^ 



For the three-point function in the region of small Gram determinants, the *4.j coefficients 



are: 



Ai 



mf — m 



2 ' 



A 2 = -A x 



1 



(38) 



We also introduce the following variables to facilitate compact expressions: 

Xij 



Pj ~ Pi 

1 2 ' 

mf — mf 
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Pj ~Pi 

2 2 

mf — mf, 

1 j 
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The zeroth order 3 J r j,(o) can be obtained from the summation of a geometric series with 
the convergence condition |%i 2 | = | X21 1 < 1: 

3 .7*1,(0) =Ai[l- X12 + X12 - X12 + ■ • ■ ] = i Al = ^ » 

1 + X21 

where Xij — Xji an d A\ = —A 2 are used. The first order 3 J r i i (i) is determined by using the 
derivative of the geometric series, 



3 J- 2 , (1) = -2a; 2 (^-^) . (40) 
The ra-th order derivative of the geometric series provides the n-th order 3 J : i 1 ( n y. 

3 T Un) = 2 n u n+i ( Pl -p 2 r^--( Pl -p 2 r- , (4i) 

3 ^2,(n) = -2 n ^ +1 ( Pl -p 2 r ■■■{p x -p 2 )^ ■ (42) 
Up to third order, Co can be compactly expressed for the case of small kinematic variables 



as: 



C = uo{Bl - B 2 ) + 2u 2 [(p? - p v p 2 )Bl + (pi - Pl p 2 )B 2 ] 
+ 4u 3 [( Pl - p 2 ) 2 {Bl - B 2 ) + {p\ - pvP2) 2 B 1 u - (p 2 2 - PiP 2 ) 2 B 2 n ] 
+ 8u 4 [3( Pl - p 2 ) 2 (pl - Pi-p 2 )Bl m + {p\ - px-p 2 fB l ni 

+ ^Pi-P2)\pl-PvP2)Bl m + {p 2 2 -p l -p 2 fB 2 lll ]+0{{p 2 l Y , (43) 

where B 1 .. is defined in Eq. (|35|) . Similarly, vector and 2 nd -rank-tensor coefficients functions 
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yield: 



d = uB\ + 2a; 2 [Bin - B 2 00 + (p 2 - prp^] 
+ 4w 3 [(3p 2 - A Pl -p 2 + pl)B] m + (p 2 - Pl ■p 2 ) 2 B\ ll 4 

+ 8c/ [3(px -P2) 2 (5 1 ooo " #0000) + 6(p?-Pi -ft) 2 ^ - 
+ 0({p\)\{p\)\{p l .p 2 f) , 



2(pl 



-prP 2 )B 2 01 ] 

(44) 



2w 2 [5, 
-4p r p 2 



2 

00 



C 2 = -wSf - 
- 4c 3 [(3p 2 

+ 8c/[3(p 1 -p 2 ) 2 (B 2 000 
+ 0{{pl)\{pl)\{p r p 2 Y) 



(pl-pi-P2)Bf 1 ] 
Pl)Bl m + (p 2 - PvP2) 2 B 2 lll + 2(p 2 - p v p 2 )Bl m ] 

) + 6(p 2 -p! -pa) 2 ^ -6(p 2 -p! -p^B 1 ^ + (pi - A -p 2 ) 3 B 2 m ] 



-^oooo 



Coo = "(Boo - 5 o 2 o) + 2^ 2 [(p 2 - PrP 2 )5 1 01 + (p 2 - p r p 2 )S 1 01 ] 
+ 4w 3 [(p! -p 2 ) 2 J B 1 000 + (p 2 -prP2) 2 ^ 011 - (pi -p 2 ) 2 B 2 000 - (p 2 -Pi-p 2 ) 2 5 2 011 ] 
+ 8w 4 [3(pi -p 2 ) 2 (p 2 -pi-p 2 )^ooooi + (Pi -PrrfSoom 

+ 3(pi - p 2 ) 2 (pi - Pi-p 2 )5 2 oooi + (P 2 - Pi-P2) 3 5o 2 om] + O ((p 2 ) 4 , (p 2 ) 4 , (p r p 2 ) 4 ) , 

(45) 

C n = o;( J B 1 1 1 ) + 2a; 2 ^Bo^ + ^-pi-^Bjn] +4c/[2B \ )00 + {4(p 2 -p 1 -p 2 ) + (p 1 -p 2 ) 2 }B \ ) 
+ (p 2 -Pi-p 2 ) 2 B 1 1 111 - 2B 2 000 ] + 8u; 4 [{6(p 2 -p r p 2 ) + 6( Pl - p 2 ) 2 }B 1 0001 
+ {6(p 2 -Pi-p 2 ) 2 + 3(pi -p 2 ) 2 (p 2 - PrP 2 ) }Boom + (p 2 -prp 2 ) 3 B 1 1 1111 
+ 6(p 2 - Pi-P 2 )B 2 0001 ] + ((p 2 ) 4 , (p 2 ) 4 , (p r p 2 ) 4 ) , 



'0011 



C12 = C 21 = -2u 2 [B^ 01 + B 2 01 ] - 8w 3 [B^ooo - B 2 000 + (pj - p 1 -p 2 )B 1 011 - (p 2 - p r p 2 )B 2 011 ] 

- 8w 4 [{6(p 2 -pi-p 2 ) + 3(pi -p 2 ) 2 }B 1 0001 + 3(p 2 -pi-p 2 ) 2 B 1 0111 

+ {6(p§ -pi-ft) + 3(p! -P 2 ) 2 }B 2 0001 + 3(pi -p 1 -p 2 ) 2 B 2 0111 ] + ((p 2 ) 4 , (p 2 ) 4 , (p r p 2 ) 4 ) , 

(46) 

C 22 =^j{B\ x )^ 2 [2B 2 01 + (p 2 - Prp 2 )B 2 n ] +4w 3 [2B 1 000 -{4(p 2 -p r p 2 ) + ( Pl -p 2 ) 2 }B 2 

- (P 2 - Pi-p 2 ) 2 B 2 lu - 2B 2 000 ] + 8a; 4 [{6(p 2 - p r p 2 ) + 6(pi - p 2 ) 2 }B 2 



•Wl 

^00001 



+ {Q{pI -Pi"P 2 ) 2 + 3(pi -p 2 ) 2 (p 2 -Pi-P 2 )}B 2 0111 + {p 2 2 - Pl -p 2 fB 2 



urn 



+ 6(p 2 -pi-p 2 )B 1 0001 ] +0((p 2 ) 4 ,(p 2 ) 4 ,(PrP 2 ) 4 ) • (47) 
Even though four-point functions are more complicated, we can apply the same technique: 
1 1 1 

/ 9 9\ / 9 9\ J *^- 2 / 9 9\ / 9 9\ J *^-3 



(raf — m 2 ){m\ — m 2 ) ' (m 2 — m\){m 2 — m 2 ) ' (m 2 — mf — m|) 



2 \ ' 
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1 



(1 + Xi 2 )(l + Xi 3 ) 



& M = - i ^ ^4 ^ va* - ( ^ ^ ^ ) r?« 



Ul1 Wl + Xl2 )2(l + Xl3 ); 12 ^ (1 + Xi2)(1+Xi3)2 ; 18. 

2 2 A \^„^ , ( 2Ai 



" + xi 2 ) 3 (i + xi,)) + V(i + xi 2 ) 2 (i + xi,?) Ws ' 

2-Ai . f 2 2 A 



(48) 



1 :i+xi 2 ) 2 (i+xi 3 ) 2 J PlW2 + l(i+xi 2 )(i+xi 3 ) 3 ]n * n * ' 



The other J-"'s can be obtained from the permutation of 1, 2 and 3 as follows: 

ijr 4tt I 4tt 4tt I 

-f~2,(n) — ■ r l,(n)|i_ >2 , 2->l ' ^3,(71)— ■ r l,(n)|i_ > .3 ) 2->l, 3->2 ' 

Explicit expressions for four-point integrals are shown in Appendix [B] Note that \Xij\ < 1 
must be satisfied. 

The numerical computations of the C and D functions are shown in both Fig. [5] and 
Fig. [7J Although Co and Dq don't have any numerical problem, the approximation of small 
kinematic variables is in accord with the analytical calculation. 

IV. SIX AND MORE POINT FUNCTIONS, F,G,H....: ALTERNATIVE REDUC- 
TION 

Here, we suggest an alternative reduction formula for the n-point integrals that are always 
decomposed to six (n — l)-point functions. Without loss of generality, pi, p 2 , P3 and p 4 
are supposed to be linearly independent. A vector pi, where t is one of the integers of 
5 < I < n — 1, can be expressed as: 

4 

pt = ^2 XiPi - 

i=l 

The be found using the 4x4 Gram matrix, Zij = 2pi-pj with i,j = 1, 2, 3, 4: 

l 



2 £(^%»'Pi • ( 49 ) 



Therefore, a system of equations, Eq. (JSJ), Eq. (J5]),and Eq. ffTU]) can be simplified as: 

a\ + aeXi = , a 2 + a.£X 2 = , 
a 3 + a^X3 = , a 4 + 0^X4 = , 
n-1 

ao + d\ H h a n _i = , ^J(p 2 - m 2 )a; = 1 . 

i=0 
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FIG. 7: Numerical safety of four-point functions : D...\pf,p^,p^, (pi — P2) 2 ,(P2 — P3) 2 ,(pi — 
p 3 ) 2 ,ml,ml,ml,ml]=D...[0,u,0,0,t,ml,M^,0,M^] with m d = 0.066 GeV. The Mandelstam 
variable, t varies from —0.0001 to 0, and s = 1 GeV 2 . The conventional Passarino-Veltman reduc- 
tion shows numerical instability in small kinematic region (red curve), while blue curve which is 
generated from the method in Sec. IIIIf Eq. (|Bip ) is numerically stable. The Feynman diagram is 
shown in Fig. [6j Red curves are generated using Loop Tools 0] . 



Hence, all coefficients a« are easily determined as: 

Xi , 



at 



U - £j fj x . 
i 



ft - Ej fo x . 



1,2,3,4) 



1 - x x - x 2 - x 3 - x 4 
ao = > ( 50 ) 

ft ~ Ef fi X 3 
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0.0 0.2 0.4 0.6 0.8 1.0 

x (=s'/s) 

FIG. 8: Comparison of six-point scalar integrals using Eq. (I14[) (or Ref. [l|) and Eq. ()5ip . The 
two curves are identical. This not only shows the validity of (]5ip but also improves numerical 
efficiency. This computation is based on the Feynman diagram in Fig. [H 

where fe=pj — m\ + ml, and I is one of the integers of 5 < I < n — 1. With the above a«'s, 
the identity for the reduction is: 

4 

1 _ aeNj y- ajNj 

jVo---jV„-i ~ ^o-JV„_i ^iVo---JVn-i ' 

Finally, an alternative reduction formula is: 

4 

T™ = a/T™- 1 ^) + OiT"- 1 ®, (n>Q), (51) 

t=0 

where I is one of the integers {5, • • • , n — 1}. The coefficients Oj are determined from the 
Xi which can be obtained from the Gram matrix (Eq. (|49p and Eq. (|50|) ). Although the 
alternative reduction of Eq. ( 15 ip uses the Gram matrix, there is no problem of the vanishing 
Gram determinant because p±, P2, P3 and p± are linearly independent. Also, the case of a 
small Gram determinant is already derived in Sec. IIHI 

According to the reduction formula of Eq. ( |T4|) . we should compute the inverse of (n + 
1) x (n + 1) Caley matrix. Here, the alternative reduction doesn't need any computation 
of a high-dimensional Caley matrix regardless of the value of n. In the previous method, 
n-point integral reduces to a sum of n number of (n — l)-point integrals. However, the 
alternative way gives always 6 number of (n — l)-point integrals which improves the efficiency 
of numerical computation. Vector and tensor coefficient functions can be obtained using the 
same technique in Sec. [TTJ 
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0.00008 
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FIG. 9: Comparison between the method in Sec lIIII (Eq. (J43]), (I46h ) and the method described in 
Ref. []J|. We computed -RefC^i^] up to second order using Eq. (H3l) . (Hoi) (dashed line), while for 
the method in Ref. three iterations used (solid curve), in which each iteration corresponds to 
the order of momentum squared. In order for the method in Ref. [l| to show reasonable stability, 
we should perform much more iterations. 



V. CONCLUSION 



We formulate a reduction for n-point integrals with a new formalism. In the region of 
small kinematic variables, which yield a small Gram determinant, it is well known that 
NLO calculation suffers from numerical instability when we use the conventional Passarino- 
Veltman reduction. The analytical expression for that region not only improves the numer- 
ical efficiency, but also fixes the numerical instability so long as the condition < 1 is 
satisfied. In Fig. and Fig. [7J we showed the numerical stability of three-point and four- 
point integrals. The numerical method described in Ref. [IJ fixes the numerical instability 
as well. In which each order of a small Gram determinant is computed iteration by iteration. 
However, the efficiency of the numerical computation is much improved as shown in Fig. [9] 
by the analytical method which provides faster convergence of the series. 

The reduction formula for the one-loop n-point integrals with n > 6, are derived alterna- 
tively in Sec. II III (Eq. floTj) ). where we reduced the number of sub-integrals. As an example, 
if we use Eq. f[T4"|) or the method suggested in Ref. [l|, then the eight-point scalar integral 
is decomposed in 8 seven-point scalar integrals. However, with the newly derived reduction 
formula, Eq. flBTl) . we needs only 6 seven-point scalar integrals. Additionally, because we 
don't use the Caley matrix, we can reduce the uncertainty which enters when computing 
the inverse of a high-dimensional Cayley matrix, thereby saving CPU time. Also, we can 
avoid the case of the vanishing Caley determinant. 
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Appendix A: Two-point integrals: B-functions 

For the 5-function expansion, we show the explicit expressions of the 5-functions and 
their derivative in this Appendix. We begin by defining the following integrals: 



J -'mn 

(a; a, /3; x , X\) 
S%, n (a;a,/3;x ,xi) 
where L 00 (a; a, (3) is: 



1 log (a(x — a)(x — P)) 

o (x — Xo) m (x — Xi) n 

1 a(x - a) m (x-(3) n 

o (x - x Q ) k (x - XiY 



(Al) 



(A2) 



L 00 (a; a, (5) = / dxlog (a(x — a)(x — (3)) = 
Jo 

( -2 - a log(l - a) + a log(-a) - (3 log(l - 0) 

+ log (a(l - a)(l - /9)) + /9 log(-/3), a, /3 ^ 0, 1 

-2-/3 log(l - /3) + /3 log(-/3) + log (a(l -/?)), a = & /3 ^ 0, 1 

-2 + log(o), « = 0&/3 = 

-2 + log(-a), « = 0&/3=l 

which is directly connected to the B function; B (pl, to.q, m\) = —L 00 (a; a, (3) with a = p\. 
The a, and /3 should satisfy the following quadratic equation: 

2 2,/ 2 2 2\ i 2 n 

p-iZ + [m l — m — p^x + m — le = . 
Because S™ e n has numerical instability in high x, let us introduce a numerically safe function 



n-l 



fn{x) 



1 1 1 hx"log (l-~) , for x < 

n- 1 n-2 1 & v x> ' 



10 



7—n 



for x > 10 



7-n 



where the number 10 7 n has been chosen for the computations shown in Fig. [5] and [7J For 
the expression of the 5-functions, we need the followings with the condition of n > — 1: 



Fi 



^rrijmia" 1 ' ^' P) 



n + 1, m, m, n + 2, 



a' /3 



/J_.q ^±£. a 0) 

LJ m,m\ a m '> ' 2 '"'^V 



a m a m /3 m (n + 1) 
1 



a m a m /3 m (n + 1) 



(0 - «)F a 



n + 1, to, m, n + 2, - , - 



a' /3 



+ 2(n + l)Fi 



n + 2, to, to, n + 3. 



i i 

a' £ 
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where Appell hypergeometric function F\ is introduced, which can be expressed in terms of 
f n using the identities of Ref. |IJ. As examples: 



™,o/i. n . .. ;m (fn(a)-fnW)) 



a(a — (5) 



0, 2M; a, 0) = — + - (f n (a) + fM) ■ 
an a 

Any order of derivatives with respect to m\ can be obtained from the following recursion 
relations: 

a qn,l ( J_. n o\ _ 1 on,0 / 1 . n o\ _ c"M-l»l f 1 ■ n 2±£. „ 

ly l J m,ml (l m ' ' 2 ' — Q w»,mVo m + 1 ' ' '"*-'m+l,m+l I a m+1 ' ' 2 ' ' ' 

where = 

The definition of L 00 gives the first derivative of the two-point scalar integral: 

C 

l a(x — a) (x — p) 



diB = — I dxdi log(a(x — a)(x — (5)) = — / dx— -, —- = — SifQ; 0; a, /?) 



o 



Although the two-point vector and tensor coefficient functions are already shown in 
Ref. (lj, we briefly derive an explicit form for them: 

= / dx(-x) n [A + log/i 2 -log(a(a;-a)(x-/3))] , (A3) 
Jo 

n 

where UV-divergence term A is defined as A = — 7_e + log(47r), and 7e is Euler's 
constant. For simplicity, the constant term, A + log/i 2 will be omitted in the following. 
After integration by parts, Eq. (IA3I) simplifies as: 



b O = -^rr ( log(a(1 ~ " ^ ~ 2aS 3? l,1 (i5 °' ^ a ' • 

n 

Now, we summarize the two-point functions and their first derivatives when p\ ^ 0: 

(_-\\n+l . x 

b O = ~^TT ( log(a(1 " " ^ " 2aS ^ ( « ; °' a ' ^J' 

A) = -L 00 (a;a,/3) 

#oo = g - log(m 2 )) + (p 2 - m 2 + m 2 )^ + 2m 2 5 + m 2 + m 2 - ±p 2 ) , 

5ooi = g ( - m\{\ - \og{m\)) + (pi -m\ + m 2 )fi n + 2m 2 5 1 - |(2m 2 + Am\ - p\) 
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n 

^o = -^i 1 ?(J;0;a,/3), 

5 00 = 1 (l - log(m 2 ) - flj + (pi -m\ + m 2 )B[ + 2m 2 £ ) , 
^ooi = g ( Mm?) - S u + (p? - m 2 + m 2 )^ + 2mg^[ - ~ 

When pj = 0: 

5 = — o— — o f m o -m" 2 -ml log(mo) + ml log(m 2 )) , 

Wq — mf V / 

51 = - m 2| 2 ( ~ 3m o + 4m o m i ~ m i + 2m o lo g( m o) ~ 2m 2 (2m 2 , - m 2 ) log(m 2 ) j , 

B oo = of 2 2T f 3m o ~ 3m i ~ 2m o lc, g( m o) + 2m i lo g( m i)) > 

Bl1 = ~ 18(m 2 -m 2 ) 3 ( 6m ° l0g ( m °) ~ 6m ?M™?) 

+ (m 2 , — ml)(— llmg + 7m.Qm, 2 — 2mf + 6(m.Q — m 2 ) 2 log(m 2 ) j , 

^ooi = 70 , 2 2^ ( - H^o + 11m\m\ - 16ml - 3m% log(mjj) + 9m^ log(ml) + 3m[j log(m 

( ^[TTIq 777. ^ j V 

Blu = 48K-mf)^ ( 12m ° l0g(m ° ) " 12m ologK) 

+ (ml - ml)(-2hml + 23mom 2 - lSm^m? + 3m? + 12(m 2 - m 2 ) 3 log(m 2 )) . 



For special cases, such as mo = or mo = mi, we have to start from the expression of 
Eq. (jSU). 
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Appendix B: D-functions for the region of small kinematic variables 

For the four-point functions described in Sec. II III for the general mass case, the explicit 
expressions in small Gram determinant region are: 

Do = a 1 B 1 +a 2 B 2 +a 3 Bl + bn(plBl- Pl - P2 Bl) + b u ( P lBl- Pl - P2 Bf)+b 13 (plBl- Pl - P3 Bl) 
+ b 21 { P \B\ - Pi-PsBl) + b 22 { P \B\ - p 2 -p 3 B\) + b 23 (p 2 3 Bf - p 2 -p 3 B 3 ) 
+ c u ((pi - P2) 2 B 1 00 + {pi - PvP^B^) + cia ((pi - p 2 ) 2 B 2 00 + {p 2 2 - Pl -p 2 ) 2 B 2 u ) 

+ C13 ((pi - p 3 ) 2 5 3 + (p| - pi -p 3 ) 2 sy 

+ 2C 2 1 ((pi - Pl -P2 - Pi 'PS + P2 -P3)-Sq + (Pi - Pi -P2) (P? ~ Pi 'Ps^n) 

+ 2c 22 ((P^ -Prp 2 +P1-P3 -p 2 -p 3 )5 2 o + (Pl -Pi-P2)(pl -P2.M) 
+ 2c 23 ((Ps +PrP2 -PrPs -P2-Ps)Bo + (pj - Pi-p 3 )(p 2 3 - p 2 -p 3 )Bf 1 ) 
+ c 3 i ((Pi - Ps) 2 ^ + (pj - Pi-PzfB 1 ^) + c 32 ((p 2 - p 3 ) 2 B 2 + (p 2 - p 2 -p 3 ) 2 B 2 u ) 

+ C33 ((P2 - P 3 ) 2 5 3 + (P3 - P2 •P 3 ) 2 5 1 3 1 ) , 

Dl = a x B\ + 6n + (p 2 - Pl -p,)^) - 6 12 5 2 - b 13 B 3 00 + 6 21 (S^ + (pj - p x -p^B^) 
+ Cii((3p 2 -4p 1 -p 2 +p^)^ 01 +(p?-p 1 -p 2 ) 2 ^ 11 )+2ci 2 (p r p 2 -p2)^ 01 +2c 13 (p 1 -p3-p^) J B^ 01 
+ 2c 2 i ((3p 2 -2(pi-p 2 +pi-p 3 ) +p 2 -p 3 )Bl m + (p 2 -prp 2 )(p 2 -Pi-p 3 )-Bj u ) 
-2c 22 (p2-p 2 -p 3 ) J Bo 01 -2c 23 (p3-p 2 -p3)^ 01 +C3i ((3p 2 -4p r p 3 +p 2 ) J Bo 01 + ^ 2 -Pi-p 3 ) 2 -Bj n ) , 

D 2 = a 2 5 2 - ftu^ij + 612 (5 2 + (p 2 - P rP2)B 2 n ) + b 22 (B 2 + (pj - p 2 -p 3 )B 2 n ) - b 23 B 3 00 
-2cii(p 2 -pi-p 2 ) J Bo 01 +ci 2 ((p 2 -4pi-p 2 + 3p 2 ) J Bo 01 +(p 2 -pi-p 2 ) 2 J B 2 11 )-2c 2 i(p 2 -pi-p 3 ) J B 1 01 
+ 2c 22 ((3p 2 -2(p r p 2 +p 2 -p 3 ) +pi-p 3 )^ 01 + (p 2 2 - Pi-p 2 )(pj - p 2 -p 3 )B 2 n ) 
-2c 23 (p 2 -prP3)^ 01 +c 32 ((3p 2 -4p 2 -p 3 +p 2 ) J Bo 01 +(p2 -p 2 -p 3 ) 2 J B 1 2 11 )+2c 3 3(p 2 -P3-P3) 5 ooi > 

L> 3 = a 3 S 3 + b 13 (B 3 + (p 2 - Pl -p 3 )B 3 u ) - b 21 B l 00 - b 22 B 2 + 6 23 (B 3 00 + (p 2 - p 2 -p 3 )Bl 1 ) 
+C13 ((p 2 -4prP3 + 3p 2 )5 3 01 + (p r p3-p^) 2 5 3 11 )-2c 21 (p 2 -prp 2 )5 V + 2c 22 (p r p 2 -p 2 )^^ 
+ 2c 23 ((p r p 2 - 2(p r p 3 + p 2 -p 3 ) + Sp 2 )^ 3 ^ + (p 2 - p r p 3 )(p 2 - p 2 -p 3 )B 3 n ) 
-2c 3 i(p 2 -pi-p 3 )5 1 01 + 2c 32 (p^-p 2 ^3)^ + 033 ((p^-4p 2 -p3+3p 2 )5 3 01 + (p 2 -p 2 -p 3 ) 2 5 1 3 11 ) . 

D 00 = a 1 B^ + a 2 B 2 + a 3 B 3 + b ll (p 2 l -p r p 2 )Bl m + b 12 (pl-p 1 -p 2 )B 2 m +b l3 (pl^ 
+ 621 (Pi -PrPs) ^ooi + MP2 -p 2 -p 3 )B 2 01 + b 23 (p 2 3 -p 2 -p 3 )B 3 0l 

+ en ((p 2 -pi-p 2 ) 2 5 1 011 + (pi - p 2 ) 2 J B 1 000 ) + c 12 ((p 2 -prP2) 2 -Boon + (pi - P2) 2 ^ 000 ) 
+ cis ((p 2 - Pi-P3) 2 5oon + (Pi " Psf^oooo) 

+ 2c 2i ((p 2 - p r p 2 )(p 2 - pvp 3 )Bl Qll + (p 2 - p r p 2 - Prp 3 + p 2 -p 3 )Bl 000 ) 

+ 2C 22 ((p 2 -p r p 2 )(p 2 -p 2 -p 3 )i?o 2 011 + {Pl -P1-P2 +PX-P3 -P2- P3 )B 2 000 ) 

+ 2c 23 ((p 2 -prp 3 )(p 2 -p 2 -p 3 )^ 011 + (pj +pi-p 2 -Pi-p 3 -p 2 -p 3 )^ 000 ) 

+ C 3 1 ((P? - Pl-p 3 ) 2 -Bo011 + (Pl - Psf^OOOo) + C 32 ((P 2 - P2-P 3 ) 2 5 2 011 + (p 2 - Pa^oOOo) 

+ C33 ((p 2 - P2-p 3 ) 2 5 3 011 + (p 2 - p 3 ) 2 B 3 000 ) , 
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D u = a.Bl, + fo n (2B l 001 + (p 2 - Pl - P2 )B 1 111 ) + b 21 (2B* 01 + (p\ - Pl -p 3 )B 1 111 ) 

+ en ((5p 2 - 6p r p 2 + P^oon + 2S oooo + (Pi - Pi-P2) 2 J B 1 1 111 ) + 2c 12j B^ 000 + 2ci 3 ^ 000 
+ 2c 2 i ((5p 2 - 3(pi-p 2 +Pi-Pa) +P2-P3)-Boon + 2Sq 000 + (p 2 - pi-p 2 )(pj - PrPsj-Bmi) 
+ c 3 i ((5p 2 - 6p r p 3 + pMoii + 2^000 + (p 2 - Prp 3 ) 2 5mi) , 

D12 = £>2i = -buB^ 01 -b 12 B 2 01 -2c u ((p 2 -p 1 •p 2 )£roii+- B oooo) -2c 12 ((p 2 2 -p x •p 2 )£roii+- B oooo) 
- 2c 2 i ((p 2 - Prp 3 )£ooii + 5 oooo) ~ 2c 22 ((P2 - ^'Ps^ooii + B% 000 ) + 2c 23 ^ 000 , 

D 13 = -613^001 -^1^001 -2ci3 {{pI -Pi -P3)£ooii + 5 oooo) -2c 2 i (2(p 2 -pi ^^oon + ^ 000 ) 
+ 2022^0000 - 2c 23 ((pg - p 2 -p 3 )Bl 0ll + ^ooo) - 2c 3 i (2(p 2 - pi-p 3 )B* 011 + 5 oooo) > 

D22 = a^n + 6i 2 (25 2 01 + (p 2 - p r p 2 )5 2 n ) + 622 (25 2 01 + (p 2 - p 2 -p 3 )B 2 ul ) + 2c 11 5 1 000 
+ C12 ((5p^ - 6pi-p 2 + p 2 )5 2 011 + 2^000 + (p 2 - Pl -p 2 ) 2 B 2 ul ) 

+ 2c 22 ((5p 2 - 3(pi-p 2 +p 2 -p 3 ) + pi-p 3 )B 2 on + 250000 + (pI - V\-Pi){p\ ~ P2-p 3 )B 2 lu ) 

+ C 3 2 ((5p 2 - 6p 2 -p 3 + P3)5 2 on + 2^000 + (P2 - P2-P 3 ) 2 5 2 m ) + 2033^000 > 

^23 = -D32 = -b 22 B 2 01 - 623^001 + 2021^0000 - 2c 22 ((pi - Pi-p 2 )B 2 011 + .Bqooo) 

- 2c 23 ((p 2 - p r p 3 )^ 011 + B 3 000 ) - 2c 32 {(p 2 2 ~ p 2 -p 3 )B 2 0U + B 2 000 ) 

- 2c 33 ((p\ - p 2 -p 3 )B 3 011 + 5 3 000 ) , 

D 33 = a 3 B 3 u + 613 (25 3 01 + (p 2 - Pl -p 3 )B 3 u ) + 6 23 (2£? 3 i + (p 2 - prp^B 3 ^) 
+ C13 {(5pl - 6prp 3 + pl)B 3 011 + 2fi 3 ooo + (Pi - Pi-P3) 2 B! 1U ) 

+ 2c 23 ((5p 2 - 3(pi-p 3 +p 2 -p 3 ) +pi-p 2 )^ 01 i + 25^00 + (pI - Pi-Ps)(pI ~ P2-p 3 )B 3 ul ) 
+ 2031^000 + 2c 32J B 2 000 + c 33 ((5p 2 - 6p 2 -p 3 + P^oon + 2b3 0Q00 + (p 2 - p 2 -p 3 ) 2 B 3 uu ) , 

(Bl) 
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with, 
ai = 
611 = 
bi2 = 
613 = 

Cll = 

C31 = 
C22 = 

Cl3 = 
C33 = 



(1 + Xi2)(l + Xis) 
1A X 



, 0,2 



A 2 



«3 



(1 + X2l)(l + X23)' (1 + X3l)(l + X32) 



(i + xi 2 ) 2 (i + xi 3 )K- 
2A 2 


m|)' 


(1 + X 2i) 2 (l + X2 3 )(m|- 
2^3 


m\) ' 


(1 + X3i) 2 (l + X3 2 )(m|- 
4A 


m 2 ) ' 


(l + Xi 2 ) 3 (l + Xi 3 )K- 


m 2 .) 2 ' 


(1 + Xi 2 )(l + Xi3) 3 (m 2 - 

2^2 


m|) 2 ' 


(i + X2i) 2 (i + x 23 ) 2 K- 

4^3 


-m 2 )( 


(1 + X 3l) 3 (l + X32)K- 
4^3 


mf) 2 ' 



^21 



'22 — 
, ^23 = 



C 2 1 = 



Cl2 



(l + Xi 2 )(l + Xi3) 2 K-m 2 )' 

2^ 

(l + X2i)(l + X23) 2 (m 2 -m 2 )' 

2^3 

4 + Xsi)(l + Xwfirni ~ ml) ' 

2^41 

(1 + Xi2) 2 (l + Xi3) 2 {ml - ml)(m\ - m§) ' 

4^ 

(l + X2i)(l + X2i) 3 (m 2 -m 2 ) 2 

4^2 

C32 



,2^ ' 



(l + X2i)(l + X2 3 ) 3 K-m 2 ) 2 ' 

2^3 

'1 + Xsi) 2 (l + X32) 2 (^l - m\){ml - ml) 



(l + X3i)(l + X3 2 ) 3 (m 2 -m 2 ) 2 



[1] A. Denner and S. Dittmaier, Nucl. Phys. B734, 62 (2006), hep-ph/0509141. 

[2] A. I. Davydychev, Phys. Lett. B263, 107 (1991). 

[3] A. Denner, Fortschr. Phys. 41, 307 (1993), 0709.1075. 

[4] Z. Bern, L. J. Dixon, and D. A. Kosower, Nucl. Phys. B412, 751 (1994), hep-ph/9306240. 

[5] T. Binoth, J. P. Guillet, and G. Heinrich, Nucl. Phys. B572, 361 (2000), hep-ph/9911342. 

[6] A. Denner and S. Dittmaier, Nucl. Phys. B658, 175 (2003), hep-ph/0212259. 

[7] G. Passarino and M. J. G. Veltman, Nucl. Phys. B160, 151 (1979). 

[8] J. Murley and N. Saad, math-ph/08095203, (2008). 

[9] T. Hahn, Nucl. Phys. Suppl. 89, 231 (2000). 

[10] G. 'tHooft and M. J. G. Veltman, Nucl. Phys. B 153, 365 (1979). 

[11] D. Binosi, J. Collins, C. Kaufhold, and L. Theussl, Comput. Phys. Commun. 180, 1709 (2009). 



24 



